Garder uniquement les OTUs présents dans >30% des échantillons (on ne fera rien des autres pour le réseau).
compter nbre de 0 dans chaque colonne et ne retenir que ceux présents dans au moins 30% des échantillons soit 253 sur 846
load("Data_Biom.Rdata")
load("MultiDataSet.Rdata")
otu_trsf <- as.data.frame(otu_table(filter_biom))
count <- apply(otu_trsf, 2, function(x){length(which(x!=0))} )
nb_count <- which(count>=253)
otu_trsf2 <- otu_trsf[,nb_count]
#as.matrix(as.integer(
#otu_table(filter_biom, taxa_are_rows = FALSE) <- otu_trsf2
otu <- as.matrix(otu_trsf2)
taxa <- tax_table(filter_biom)
Sample_ID <- sample_data(filter_biom)
Fil_Biom <- phyloseq(otu_table(otu, taxa_are_rows = FALSE),
tax_table(taxa),
sample_data(Sample_ID))
taxa <- as.data.frame(tax_table(Fil_Biom))
filter_biom <- Fil_BiomSortie : >huge::huge.roc(se\(est\)path, graph, verbose=FALSE)
True Postive Rate: from 0 to 0.988806 False Positive Rate: from 0 to 0.7952971 Area under Curve: 0.7800105 Maximum F1 Score: 0.9760311 > stars.pr(getOptMerge(se), graph, verbose=FALSE)
True Postive Rate: from 0 to 0.9626866 False Positive Rate: from 0 to 0.1032948 Area under Curve: 0.9778643 Maximum F1 Score: 0.9642207
MB, Glasso et SparCC
#load("spiec_easi.Rdata")
#load("spiec_easi2.Rdata")
#load("spiec_easi_filt.Rdata")
#se.mb.biom_init <- se.mb.biom
#se.gl.biom_init <- se.gl.biom
#load("spiec_easi_filt_teste1.Rdata")
#load("spiec_easi_filt_teste2.Rdata")
#load("spiec_easi_p1.Rdata")
#load("spiec_easi_gl_mb.Rdata")
load("spiec_easi_gl_mb.Rdata")
## Define arbitrary threshold for SparCC correlation matrix for the graph
sparcc.graph <- abs(sparcc.biom$Cor) >= 0.3
diag(sparcc.graph) <- 0
sparcc.graph <- Matrix(sparcc.graph, sparse=TRUE)
## Create igraph objects
#lambda.min.ratio=1e-3
ig.mb_1 <- adj2igraph(getRefit(se.mb.biom_1))
ig.gl_1 <- adj2igraph(getRefit(se.gl.biom_1))
ig.sparcc_1 <- adj2igraph(sparcc.graph)
#lambda.min.ratio=1e-2 val par défaut
ig.mb <- adj2igraph(getRefit(se.mb.biom))
ig.gl <- adj2igraph(getRefit(se.gl.biom))
ig.sparcc <- adj2igraph(sparcc.graph)
#lambda.min.ratio=5e-2
ig.mb005 <- adj2igraph(getRefit(se.mb.biom005))
ig.gl005 <- adj2igraph(getRefit(se.gl.biom005))
ig.sparcc005 <- adj2igraph(sparcc.graph)
#lambda.min.ratio=0.1
ig.mb01 <- adj2igraph(getRefit(se.mb.biom01))
ig.gl01 <- adj2igraph(getRefit(se.gl.biom01))
ig.sparcc01 <- adj2igraph(sparcc.graph)
#lambda.min.ratio=0.15
ig.mb015 <- adj2igraph(getRefit(se.mb.biom015))
ig.gl015 <- adj2igraph(getRefit(se.gl.biom015))
ig.sparcc015 <- adj2igraph(sparcc.graph)
#lambda.min.ratio=0.2
ig.mb02 <- adj2igraph(getRefit(se.mb.biom02))
ig.gl02 <- adj2igraph(getRefit(se.gl.biom02))
ig.sparcc02 <- adj2igraph(sparcc.graph)
#lambda.min.ratio=0.3
ig.mb03 <- adj2igraph(getRefit(se.mb.biom03))
ig.gl03 <- adj2igraph(getRefit(se.gl.biom03))
ig.sparcc03 <- adj2igraph(sparcc.graph)## set size of vertex proportional to clr-mean
vsize <- rowMeans(clr(otu_table(filter_biom), 1))+6
biom.coord_1 <- layout.fruchterman.reingold(ig.mb_1)
biom.coord <- layout.fruchterman.reingold(ig.mb)
biom.coord005 <- layout.fruchterman.reingold(ig.mb005)
biom.coord01 <- layout.fruchterman.reingold(ig.mb01)
biom.coord015 <- layout.fruchterman.reingold(ig.mb015)
biom.coord02 <- layout.fruchterman.reingold(ig.mb02)
biom.coord03 <- layout.fruchterman.reingold(ig.mb03)
biom.coord_g1 <- layout.fruchterman.reingold(ig.gl_1)
biom.coordg <- layout.fruchterman.reingold(ig.gl)
biom.coordg005 <- layout.fruchterman.reingold(ig.gl005)
biom.coordg01 <- layout.fruchterman.reingold(ig.gl01)
biom.coordg015 <- layout.fruchterman.reingold(ig.gl015)
biom.coordg02 <- layout.fruchterman.reingold(ig.gl02)
biom.coordg03 <- layout.fruchterman.reingold(ig.gl03)
sparcc.graph1 <- abs(sparcc.biom$Cor) >= 0.01
diag(sparcc.graph) <- 0
sparcc.graph1 <- Matrix(sparcc.graph, sparse=TRUE)
ig.sparcct1 <- adj2igraph(sparcc.graph1)
sparcc.graph2 <- abs(sparcc.biom$Cor) >= 0.9
diag(sparcc.graph) <- 0
sparcc.graph2 <- Matrix(sparcc.graph, sparse=TRUE)
ig.sparcct2 <- adj2igraph(sparcc.graph2)
par(mfrow=c(1,3))
plot(ig.sparcc, layout=biom.coord, vertex.size=vsize, vertex.label=NA, main="sparcc 0.3 layoutmb")
plot(ig.sparcct1, layout=biom.coord, vertex.size=vsize, vertex.label=NA, main="sparcc 0.01")
plot(ig.sparcct2, layout=biom.coord, vertex.size=vsize, vertex.label=NA, main="sparcc 0.9")par(mfrow=c(1,3))
plot(ig.sparcc, layout=biom.coordg, vertex.size=vsize, vertex.label=NA, main="sparcc 0.3 layoutglasso")
plot(ig.sparcct1, layout=biom.coordg, vertex.size=vsize, vertex.label=NA, main="sparcc 0.01")
plot(ig.sparcct2, layout=biom.coordg, vertex.size=vsize, vertex.label=NA, main="sparcc 0.9")par(mfrow=c(1,3))
plot(ig.mb_1, layout=biom.coord_1, vertex.size=vsize, vertex.label=NA, main="MB λ=1e-3 ")
plot(ig.gl_1, layout=biom.coord_g1, vertex.size=vsize, vertex.label=NA, main="glasso")
plot(ig.sparcc_1, layout=biom.coord_1, vertex.size=vsize, vertex.label=NA, main="sparcc")par(mfrow=c(1,3))
plot(ig.mb, layout=biom.coord, vertex.size=vsize, vertex.label=NA, main="MB λ=1e-2 ")
plot(ig.gl, layout=biom.coordg, vertex.size=vsize, vertex.label=NA, main="glasso")
plot(ig.sparcc, layout=biom.coord, vertex.size=vsize, vertex.label=NA, main="sparcc")par(mfrow=c(1,3))
plot(ig.mb005, layout=biom.coord005, vertex.size=vsize, vertex.label=NA, main="MB λ=0.05 ")
plot(ig.gl005, layout=biom.coordg005, vertex.size=vsize, vertex.label=NA, main="glasso")
plot(ig.sparcc005, layout=biom.coord005, vertex.size=vsize, vertex.label=NA, main="sparcc")par(mfrow=c(1,3))
plot(ig.mb01, layout=biom.coord01, vertex.size=vsize, vertex.label=NA, main="MB λ=0.1 ")
plot(ig.gl01, layout=biom.coordg01, vertex.size=vsize, vertex.label=NA, main="glasso")
plot(ig.sparcc01, layout=biom.coord01, vertex.size=vsize, vertex.label=NA, main="sparcc")par(mfrow=c(1,3))
plot(ig.mb015, layout=biom.coord015, vertex.size=vsize, vertex.label=NA, main="MB λ=0.15 ")
plot(ig.gl015, layout=biom.coordg015, vertex.size=vsize, vertex.label=NA, main="glasso")
plot(ig.sparcc, layout=biom.coord015, vertex.size=vsize, vertex.label=NA, main="sparcc")par(mfrow=c(1,3))
plot(ig.mb02, layout=biom.coord02, vertex.size=vsize, vertex.label=NA, main="MB λ=0.2 ")
plot(ig.gl02, layout=biom.coordg02, vertex.size=vsize, vertex.label=NA, main="glasso")
plot(ig.sparcc02, layout=biom.coord02, vertex.size=vsize, vertex.label=NA, main="sparcc")par(mfrow=c(1,3))
plot(ig.mb03, layout=biom.coord03, vertex.size=vsize, vertex.label=NA, main="MB λ=0.3 ")
plot(ig.gl03, layout=biom.coordg03, vertex.size=vsize, vertex.label=NA, main="glasso")
plot(ig.sparcc03, layout=biom.coord03, vertex.size=vsize, vertex.label=NA, main="sparcc")secor_1 <- cov2cor(getOptCov(se.gl.biom_1))
sebeta_1 <- symBeta(getOptBeta(se.mb.biom_1), mode='maxabs')
elist.gl_1 <- summary(triu(secor_1*getRefit(se.gl.biom_1), k=1))
elist.mb_1 <- summary(sebeta_1)
elist.sparcc <- summary(sparcc.graph*sparcc.biom$Cor)
hist(elist.mb_1[,3], main='λ=1e-3 grey=MB\n green=sparcc, red=Glasso', xlab='edge weights')
hist(elist.sparcc[,3], add=TRUE, col='forestgreen')
hist(elist.gl_1[,3], add=TRUE, col='red')secor <- cov2cor(getOptCov(se.gl.biom))
sebeta <- symBeta(getOptBeta(se.mb.biom), mode='maxabs')
elist.gl <- summary(triu(secor*getRefit(se.gl.biom), k=1))
elist.mb <- summary(sebeta)
elist.sparcc <- summary(sparcc.graph*sparcc.biom$Cor)
hist(elist.mb[,3], main='λ=1e-2 grey=MB\n green=sparcc, red=Glasso', xlab='edge weights')
hist(elist.sparcc[,3], add=TRUE, col='forestgreen')
hist(elist.gl[,3], add=TRUE, col='red')secor005 <- cov2cor(getOptCov(se.gl.biom005))
sebeta005 <- symBeta(getOptBeta(se.mb.biom005), mode='maxabs')
elist.gl005 <- summary(triu(secor005*getRefit(se.gl.biom005), k=1))
elist.mb005 <- summary(sebeta005)
elist.sparcc <- summary(sparcc.graph*sparcc.biom$Cor)
hist(elist.mb005[,3],ylim=c(1,400), main='λ= 0.05 grey=MB\n green=sparcc, red=Glasso', xlab='edge weights')
hist(elist.sparcc[,3], add=TRUE, col='forestgreen')
hist(elist.gl005[,3], add=TRUE, col='red')secor01 <- cov2cor(getOptCov(se.gl.biom01))
sebeta01 <- symBeta(getOptBeta(se.mb.biom01), mode='maxabs')
elist.gl01 <- summary(triu(secor005*getRefit(se.gl.biom01), k=1))
elist.mb01 <- summary(sebeta01)
elist.sparcc <- summary(sparcc.graph*sparcc.biom$Cor)
hist(elist.mb01[,3], main='λ=0.1 grey=MB\n green=sparcc, red=Glasso', xlab='edge weights')
hist(elist.sparcc[,3], add=TRUE, col='forestgreen')
hist(elist.gl01[,3], add=TRUE, col='red')secor015 <- cov2cor(getOptCov(se.gl.biom015))
sebeta015 <- symBeta(getOptBeta(se.mb.biom015), mode='maxabs')
elist.gl015 <- summary(triu(secor015*getRefit(se.gl.biom015), k=1))
elist.mb015 <- summary(sebeta015)
elist.sparcc <- summary(sparcc.graph*sparcc.biom$Cor)
hist(elist.mb015[,3], main='λ=0.15 grey=MB\n green=sparcc, red=Glasso', xlab='edge weights')
hist(elist.sparcc[,3], add=TRUE, col='forestgreen')
hist(elist.gl015[,3], add=TRUE, col='red')secor02 <- cov2cor(getOptCov(se.gl.biom02))
sebeta02 <- symBeta(getOptBeta(se.mb.biom02), mode='maxabs')
elist.gl02 <- summary(triu(secor02*getRefit(se.gl.biom02), k=1))
elist.mb02 <- summary(sebeta02)
elist.sparcc <- summary(sparcc.graph*sparcc.biom$Cor)
hist(elist.mb02[,3], main='λ=0.2 grey=MB\n green=sparcc, red=Glasso', xlab='edge weights')
hist(elist.sparcc[,3], add=TRUE, col='forestgreen')
hist(elist.gl02[,3], add=TRUE, col='red')secor03 <- cov2cor(getOptCov(se.gl.biom03))
sebeta03 <- symBeta(getOptBeta(se.mb.biom03), mode='maxabs')
elist.gl03 <- summary(triu(secor*getRefit(se.gl.biom03), k=1))
elist.mb03 <- summary(sebeta03)
elist.sparcc <- summary(sparcc.graph*sparcc.biom$Cor)
hist(elist.mb03[,3], main='λ=0.3 grey=MB\n green=sparcc, red=Glasso', xlab='edge weights')
hist(elist.sparcc[,3], add=TRUE, col='forestgreen')
hist(elist.gl03[,3], add=TRUE, col='red')dd.gl01 <- degree.distribution(ig.gl01)
dd.mb01 <- degree.distribution(ig.mb01)
dd.sparcc <- degree.distribution(ig.sparcc)
plot(0:(length(dd.sparcc)-1), dd.sparcc, ylim=c(0,.35), type='b',
ylab="Frequency", xlab="Degree", main="Degree Distributions λ=0.1 ")
points(0:(length(dd.gl01)-1), dd.gl01, col="red" , type='b')
points(0:(length(dd.mb01)-1), dd.mb01, col="forestgreen", type='b')
legend("topright", c("MB", "glasso", "sparcc"),
col=c("forestgreen", "red", "black"), pch=1, lty=1)dd.gl015 <- degree.distribution(ig.gl015)
dd.mb015 <- degree.distribution(ig.mb015)
dd.sparcc <- degree.distribution(ig.sparcc)
plot(0:(length(dd.sparcc)-1), dd.sparcc, ylim=c(0,.35), type='b',
ylab="Frequency", xlab="Degree", main="Degree Distributions λ=0.15")
points(0:(length(dd.gl015)-1), dd.gl015, col="red" , type='b')
points(0:(length(dd.mb015)-1), dd.mb015, col="forestgreen", type='b')
legend("topright", c("MB", "glasso", "sparcc"),
col=c("forestgreen", "red", "black"), pch=1, lty=1)dd.gl02 <- degree.distribution(ig.gl02)
dd.mb02 <- degree.distribution(ig.mb02)
dd.sparcc <- degree.distribution(ig.sparcc)
plot(0:(length(dd.sparcc)-1), dd.sparcc, ylim=c(0,.35), type='b',
ylab="Frequency", xlab="Degree", main="Degree Distributions λ=0.2")
points(0:(length(dd.gl02)-1), dd.gl02, col="red" , type='b')
points(0:(length(dd.mb02)-1), dd.mb02, col="forestgreen", type='b')
legend("topright", c("mb", "glasso", "sparcc"),
col=c("forestgreen", "red", "black"), pch=1, lty=1)## Load round 2 of American gut project
# data('amgut2.filt.phy')
# table_ex <- as.data.frame(tax_table(amgut2.filt.phy))
# se.mb.amgut2 <- spiec.easi(filter_biom, method='mb', lambda.min.ratio=1e-2,
# nlambda=20, pulsar.params=list(rep.num=50))
# table_otu <- as.data.frame(tax_table(filter_biom))
#load("param_spiec_clut1.Rdata")
#load("param_spiec_filt2.Rdata")
#lambda.min.ratio=0.1
biom.coord_1 <- layout.fruchterman.reingold(ig.mb_1)
biom.coord <- layout.fruchterman.reingold(ig.mb)
biom.coord005 <- layout.fruchterman.reingold(ig.mb005)
biom.coord01 <- layout.fruchterman.reingold(ig.mb01)
biom.coord015 <- layout.fruchterman.reingold(ig.mb015)
biom.coord02 <- layout.fruchterman.reingold(ig.mb02)
biom.coord03 <- layout.fruchterman.reingold(ig.mb03)
ig_1.mb_1 <- adj2igraph(getRefit(se.mb.biom_1), vertex.attr=list(name=taxa_names(filter_biom)))
plot_network(ig_1.mb_1, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="λ=1e-3" )ig_1.mb <- adj2igraph(getRefit(se.mb.biom), vertex.attr=list(name=taxa_names(filter_biom)))
plot_network(ig_1.mb, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="λ=1e-2")ig_1.mb005 <- adj2igraph(getRefit(se.mb.biom005), vertex.attr=list(name=taxa_names(filter_biom)))
plot_network(ig_1.mb005, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="λ=0.005")ig_1.mb01 <- adj2igraph(getRefit(se.mb.biom01), vertex.attr=list(name=taxa_names(filter_biom)))
plot_network(ig_1.mb01, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="λ=0.1")ig_1.mb015 <- adj2igraph(getRefit(se.mb.biom015), vertex.attr=list(name=taxa_names(filter_biom)))
plot_network(ig_1.mb015, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="λ=0.15")ig_1.mb02 <- adj2igraph(getRefit(se.mb.biom02), vertex.attr=list(name=taxa_names(filter_biom)))
plot_network(ig_1.mb01, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="λ=0.2")ig_1.mb03 <- adj2igraph(getRefit(se.mb.biom03), vertex.attr=list(name=taxa_names(filter_biom)))
plot_network(ig_1.mb03, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="λ=0.3")# library(phyloseq)
# ## Load round 2 of American gut project
# data('amgut2.filt.phy')
# se.mb.amgut2 <- spiec.easi(amgut2.filt.phy, method='mb', lambda.min.ratio=1e-2,
# nlambda=20, pulsar.params=list(rep.num=50))
# ig2.mb <- adj2igraph(getRefit(se.mb.amgut2), vertex.attr=list(name=taxa_names(amgut2.filt.phy)))
# plot_network(ig2.mb, amgut2.filt.phy, type='taxa', color="Rank3")biom.coord_1 <- layout.fruchterman.reingold(ig.gl_1)
biom.coord <- layout.fruchterman.reingold(ig.gl)
biom.coord005 <- layout.fruchterman.reingold(ig.gl005)
biom.coord01 <- layout.fruchterman.reingold(ig.gl01)
biom.coord015 <- layout.fruchterman.reingold(ig.gl015)
biom.coord02 <- layout.fruchterman.reingold(ig.gl02)
biom.coord03 <- layout.fruchterman.reingold(ig.gl03)
ig_1.gl_1 <- adj2igraph(getRefit(se.gl.biom_1), vertex.attr=list(name=taxa_names(filter_biom)))
plot_network(ig_1.gl_1, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="λ=1e-3" )ig_1.gl <- adj2igraph(getRefit(se.gl.biom), vertex.attr=list(name=taxa_names(filter_biom)))
plot_network(ig_1.gl, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="λ=1e-2")ig_1.gl005 <- adj2igraph(getRefit(se.gl.biom005), vertex.attr=list(name=taxa_names(filter_biom)))
plot_network(ig_1.gl005, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="λ=0.005")ig_1.gl01 <- adj2igraph(getRefit(se.gl.biom01), vertex.attr=list(name=taxa_names(filter_biom)))
plot_network(ig_1.gl01, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="λ=0.1")ig_1.gl015 <- adj2igraph(getRefit(se.gl.biom015), vertex.attr=list(name=taxa_names(filter_biom)))
plot_network(ig_1.gl015, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="λ=0.15")ig_1.gl02 <- adj2igraph(getRefit(se.gl.biom02), vertex.attr=list(name=taxa_names(filter_biom)))
plot_network(ig_1.gl01, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="λ=0.2")ig_1.gl03 <- adj2igraph(getRefit(se.gl.biom03), vertex.attr=list(name=taxa_names(filter_biom)))
plot_network(ig_1.gl03, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="λ=0.3")load("param_spiec_slr.Rdata")
#se.slr10.biom <- spiec.easi(otu_table(filter_biom), method='slr', r=10, lambda.min.ratio=1e-2,
# nlambda=20, pulsar.params=list(rep.num=50, ncores=16))
#Différents param testés slr 10 et 20 lambda e-2 et 0.1
ig.slr10 <- adj2igraph(getRefit(se.slr10.biom))
biom.coordslr10 <- layout.fruchterman.reingold(ig.slr10)
ig.slr20 <- adj2igraph(getRefit(se1.slr20.biom))
biom.coordslr20 <- layout.fruchterman.reingold(ig.slr20)
ig.slr101 <- adj2igraph(getRefit(se.slr10.biom1))
biom.coordslr101 <- layout.fruchterman.reingold(ig.slr101)
ig.slr201 <- adj2igraph(getRefit(se1.slr20.biom1))
biom.coordslr201 <- layout.fruchterman.reingold(ig.slr201)
vsize <- rowMeans(clr(otu_table(filter_biom), 1))+6
par(mfrow=c(2,2))
plot(ig.slr10, layout=biom.coordslr10, vertex.size=vsize, vertex.label=NA, main="slr 10 λ=1e-2")
plot(ig.slr20, layout=biom.coordslr20, vertex.size=vsize, vertex.label=NA, main="slr 20 λ=1e-2")
plot(ig.slr101, layout=biom.coordslr101, vertex.size=vsize, vertex.label=NA, main="slr 10 λ=0.1")
plot(ig.slr101, layout=biom.coordslr201, vertex.size=vsize, vertex.label=NA, main="slr 20 λ=0.1")ig1.slr10 <- adj2igraph(getRefit(se.slr10.biom), vertex.attr=list(name=taxa_names(filter_biom)))
plot_network(ig1.slr10, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="slr 10 λ=1e-2" )ig1.slr20 <- adj2igraph(getRefit(se1.slr20.biom), vertex.attr=list(name=taxa_names(filter_biom)))
plot_network(ig1.slr20, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="slr 20 λ=1e-2" )ig1.slr101 <- adj2igraph(getRefit(se.slr10.biom1), vertex.attr=list(name=taxa_names(filter_biom)))
plot_network(ig1.slr101, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="slr 10 λ=0.1" )